// ****************************************************************************
// *     Filename:  PrimaryGeneratorAction.cc
// *      Created:  2012-09-14 11:58:25
// *      Version:  1.0
// *     Revision:  0
// *       Author:  Tran Hoai Nam (nam), nam@kuno-g.phys.sci.osaka-u.ac.jp
// * Organization:  Osaka University
// *      License:  Geant4 Software License 
// *                http://geant4.cern.ch/license/LICENSE.html
// *  Description:  
// ****************************************************************************

#include "PrimaryGeneratorAction.hh"
#include "DetectorConstruction.hh"

#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "Randomize.hh"
#include "globals.hh"

#include "CLHEP/Random/RandFlat.h"

#include "iostream"

PrimaryGeneratorAction::PrimaryGeneratorAction() 
{; }

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  delete fpGun;
}

PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction *det) 
	: fpDet(det)
{
  G4int n_particle = 1;
  fpGun = new G4ParticleGun(n_particle);
  fpTable = G4ParticleTable::GetParticleTable();
}

void PrimaryGeneratorAction::SetupParticleGun() 
{
	// Primary particle type
  fpGun->SetParticleDefinition(fpTable->FindParticle("mu-"));
  fpGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));

	// Primary position 
	G4double posZ = -0.5*fpDet->GetWorldZ();
	G4double posX = GaussianPos(0,5)*mm;
	G4double posY = GaussianPos(0,5)*mm;
	//G4double posX = FlatPos(-10,10)*mm;
	//G4double posY = FlatPos(-10,10)*mm;
  fpGun->SetParticlePosition(G4ThreeVector(posX,posY,posZ));

	// Primary energy
	G4double momentum = GaussianMomentum(32,2/100.);
	G4double mass = fpGun->GetParticleDefinition()->GetPDGMass()/MeV;
	G4double KE = sqrt(mass*mass + momentum*momentum) - mass;
	fpGun->SetParticleEnergy(KE*MeV);
}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
	int n_pri = 1;
	for (int i = 0; i < n_pri; i++) 
	{
		SetupParticleGun();
		fpGun->GeneratePrimaryVertex(anEvent);
	}
}

// Generate gaussianly random momentum with a centroid and relative fwhm
// Input: central momentum & fwhm
// Output: randomized (Gaussian) momentum; Unit is MeV/c
G4double PrimaryGeneratorAction::GaussianMomentum(G4double p0, G4double fwhm) 
{
	G4double sigma = fwhm*p0/2.35;
	return CLHEP::RandGauss::shoot(p0,sigma);
}
// Generate random primary position, flat distribution between xmin-xmax and
// ymin - ymax
G4double PrimaryGeneratorAction::FlatPos(G4double xmin, G4double xmax)
{
	return CLHEP::RandFlat::shoot(xmin,xmax);
}

G4double PrimaryGeneratorAction::GaussianPos(G4double x0, G4double sigma) 
{
	return CLHEP::RandGauss::shoot(x0,sigma);
}
//
G4double PrimaryGeneratorAction::GetPrimaryEnergy() 
{
	return fpGun->GetParticleEnergy();
}

G4ThreeVector PrimaryGeneratorAction::GetPrimaryPos() 
{
	return fpGun->GetParticlePosition();
}

